Data Analysis of Aviation Safety Occurance

Data Analysis Project Phase 1: Data Gathering and Exploratory Data Analysis

Amirabbas Jalali - Mina Moosavifar


در این پروژه رخدادهای امنیت پروازها از سال ۱۹۱۹ تا کنون را بررسی می کنیم.


توضیحات نحوه ی جمع آوری داده ها و پاکسازی آن در فایل report.pdf آمده است.
داده های ذخیره شده از طریق لینک زیر قابل دسترسی است:
https://github.com/Ajal88/DA_Project


کتابخانه های مورد نیاز

library(readr)
library(dplyr)
library(stringr)
library(highcharter)
library(ggplot2)
library(stringr)
library(topicmodels)
library(tidytext)

پاکسازی داده ها

توضیحات پیاده سازی پاکسازی داده ها در فایل گزارش آمده است.

asn <- readr::read_csv("Data/asn.csv")

# delete junk rows
asn <- asn %>% arrange(Engines)
asn = asn[-c(1:2),]

# fix shifted data
bug <- asn %>% arrange(desc(Engines)) %>% slice(1:10)
asn <- asn %>% arrange(desc(Engines))
asn = asn[-c(1:12),]
asn$Engines = as.numeric(asn$Engines)
names_list <- colnames(bug)
x = names_list[1:3]
y = names_list[5:21]
x = append(x,y)
x[21] = "Operator"
colnames(bug) = x
bug <- bug %>% select(colnames(asn))
bug$Engines = as.numeric(bug$Engines)
asn <- bind_rows(asn, bug)
remove(bug, x, y, names_list)

asn <- asn %>% select(-Time, -Registration, -C.n.msn)

# change strings to data for passengers
s <- asn %>% select(Passengers)
s <- str_split_fixed(s$Passengers, pattern = " ", 5) %>% as.data.frame(stringsAsFactors = F)
colnames(s) = c("kill","kill_count", "slash", "busy", "busy_kill_count")
s <- s %>% select(Passenger_fatalities = kill_count, Passenger_occupants = busy_kill_count)
s[is.na(s)] <- 0
s[s == ""] <- 0
s$Passenger_fatalities = as.numeric(s$Passenger_fatalities)
s$Passenger_occupants = as.numeric(s$Passenger_occupants)

asn <- asn %>% select(-Passengers) %>% append(s) %>% as.data.frame(stringsAsFactors = F)
remove(s)

# change strings to data for crew
s <- str_split_fixed(asn$Crew, pattern = " ", 5) %>% as.data.frame(stringsAsFactors = F)
colnames(s) = c("kill","kill_count", "slash", "busy", "busy_kill_count")
s <- s %>% select(Crew_fatalities = kill_count, Crew_occupants = busy_kill_count)
s[is.na(s)] <- 0
s[s == ""] <- 0
s$Crew_fatalities = as.numeric(s$Crew_fatalities)
s$Crew_occupants = as.numeric(s$Crew_occupants)

asn <- asn %>% select(-Crew) %>% append(s) %>% as.data.frame(stringsAsFactors = F)
remove(s)

# change strings to data for passengers
asn <- asn %>% as.data.frame(stringsAsFactors = F) %>% select(-Total) %>% 
  mutate(Total_fatalities = Crew_fatalities + Passenger_fatalities,
         Total_occupants = Crew_occupants + Passenger_occupants)

# change date of first flight
asn = asn[-c(17376),]

asn$FirstFlight <- str_remove_all(asn$FirstFlight, "\\s")

asn <- asn %>% mutate(find = grepl(pattern = "^\\d{4}", asn$FirstFlight))
asn_false <- asn %>% filter(find == FALSE) %>% mutate(FirstFlight = NA)
asn_true <- asn %>% filter(find == TRUE)
asn <- asn %>% select(-FirstFlight)
asn_true$FirstFlight = str_match(asn_true$FirstFlight, pattern = "^\\d{4}")
asn <- bind_rows(asn_false, asn_true) %>% select(-find)

remove(asn_false, asn_true)

# change date of first flight
asn <- asn %>% mutate(find = grepl(pattern = "\\d{4}", asn$Date))
asn_true <- asn %>% filter(find == TRUE)
asn_true$Date = str_match(asn_true$Date, pattern = "\\d{4}")
asn_false <- asn %>% filter(find == FALSE) %>% mutate(Date = NA)
asn <- asn %>% select(-Date)
asn <- bind_rows(asn_false, asn_true)
asn <- asn %>% select(-find)

remove(asn_false, asn_true)

# change data of aeroflot
asn <- asn %>% mutate(find = str_detect(Operator, regex("aeroflot",ignore_case = T)))
asn_true <- asn %>% filter(find == TRUE)
asn_true$Operator = "Aeroflot"
asn_false <- asn %>% filter(find == FALSE)
asn <- bind_rows(asn_false, asn_true)
asn <- asn %>% select(-find)

remove(asn_false, asn_true)

asn[asn == "Unknown"] <- NA

asn[asn == "unknown"] <- NA

asn[asn == ""] <- NA

asn <- asn %>% mutate(DepartureAirport = ifelse(DepartureAirport == "?" | DepartureAirport == "-", NA, DepartureAirport))

write.csv(asn, file = "Data/asn_c.csv",row.names=FALSE)

در نهایت داده در فایل asn_c.csv قرار دارد.


بارگذاری داده
برای سهولت کار با داده، ستون بازماندگان و نرخ زنده ماندن را به داده ها اضافه می کنیم. همچنین ایرادات این داده را که شامل خطای وارد شدن تعداد مسافران و کشته ها است را نیز برطرف می کنیم. هم چنین ستونی به عنوان نظامی بودن هواپیما و خط هوایی را نیز محاسبه می کنیم. در نهایت نیز ستون شماره ی تصادف را اضافه می کنیم که این ستون بعدا برای دسته بندی استفاده می شود.

casn <- readr::read_csv("Data/asn_c.csv") %>% as.data.frame(stringsAsFactors = F) %>% 
  mutate(Total_occupants = ifelse(Total_occupants == 0 & Total_fatalities != 0, Total_fatalities, Total_occupants),
         Total_survivors = abs(Total_occupants - Total_fatalities)) %>% 
  mutate(Total_survivors = ifelse(Total_survivors > Total_occupants, Total_occupants, abs(Total_occupants - Total_fatalities))) %>% 
  mutate(is_army = str_detect(Operator, regex("Force|Navy",ignore_case = T))) %>% 
  mutate(occ_no = row_number())

بررسی روند تلفات رخدادهای امنیتی پروازها در طول سالیان
برای این منظور ابتدا داده ها را بر اساس سال گروه بندی می کنیم، سپس تعداد تلفات، بازماندگان، افراد درگیر در حادثه و نرخ زنده ماندن را بدست می آوریم.(برای سال ۱۹۲۱ داده ی مناسبی به وجود نداشت به همین علت این سال از داده ها حذف شده است.)

# army and civil flights
year_fat <- casn %>% filter(!is.na(Date)) %>% group_by(Date) %>% 
  summarise(Total_occupants = sum(Total_occupants), Total_fatalities = sum(Total_fatalities),
            Total_survivors = sum(Total_survivors), Survival_rate = 100*Total_survivors/Total_occupants)

# remove bad data
year_fat <- year_fat[-c(3),]

highchart() %>% 
  hc_add_series(data = year_fat, type = "spline", hcaes(x = Date, y = Total_fatalities), name = "Total Fatalities") %>% 
  hc_add_series(data = year_fat, type = "spline", hcaes(x = Date, y = Total_survivors), name = "Total Survivors") %>% 
  hc_yAxis(title = list(text = "Count")) %>% 
  hc_xAxis(title = list(text = "Year")) %>% 
  hc_title(text = "Fatalities Per Year", style = list(fontWeight = "bold")) %>%
  hc_add_theme(hc_theme_flat())

همانطور که مشاهده می کنیم، تلفات حوادث در حال کاهش است. البته باید دقت داشته باشیم که این کاهش هم چنین نشانگر این است که استاندارد پرواز ها بالاتر رفته است. زیرا هر چه سال جلوتر می روند، تکنولوژی نیز پیشرفت کرده و تعداد مسافران هواپیماها افزایش یافته و استفاده از سفر هوایی بیشتر می شود. پس تعداد مسافرین بیشتر شده و تعداد کشتگان کمتر می شود که نشان دهنده ی بهبود وضعیت است.

year_fat  %>% 
  hchart(type = "spline", hcaes(x = Date, y = Survival_rate), name = "Survival Rate") %>% 
  hc_yAxis(title = list(text = "Survival Rate")) %>% 
  hc_xAxis(title = list(text = "Year")) %>% 
  hc_title(text = "Survival Rate Per Year", style = list(fontWeight = "bold")) %>%
  hc_add_theme(hc_theme_sandsignika())

با توجه به نمودارهای بالا همانطور که انتظار داشتیم، نرخ زنده ماندن تقریبا به صورت خطی بیشتر شده است.


بدترین خطوط هوایی، بدترین هواپیماها و بدترین فرودگاه ها
ابتدا معیار بد بودن را انتخاب می کنیم، از آنجایی که شرکت هایی که تعداد پایینی پرواز و یا پروازهای کوچکی داشته باشند، در صورت سقوط دارای نرخ پایین زنده ماندن هستند اما در واقع حادثه بزرگی به شمار نمی آیند، تنها خطوط هوایی، هواپیماها و فرودگاه هایی را انتخاب می کنیم که بالای ۵۰۰ نفر مسافر داشته اند. سپس معیار بد بودن را نرخ پایین زنده ماندن در نظر میگیریم.

# worst airline
worst_airline <- casn %>% filter(!is.na(Operator)) %>% 
  filter(is_army == FALSE) %>% 
  group_by(Operator) %>% 
  summarise(Total_occupants = sum(Total_occupants), Total_fatalities = sum(Total_fatalities),
            Total_survivors = sum(Total_survivors), Survival_rate = 100*Total_survivors/Total_occupants) %>% 
  ungroup() %>% 
  filter(Total_occupants > 500) %>% 
  top_n(20, wt = desc(Survival_rate)) %>% 
  arrange(Survival_rate)

p = ggplot(data = worst_airline, mapping = aes(x = reorder(Operator, Survival_rate), y = Survival_rate, fill = Total_fatalities)) + 
  geom_bar(stat="identity") + scale_fill_gradient(low="brown1", high="brown4") + 
  ggtitle("Worst Airlines with lowest survival rate") + 
  xlab("Airline") + 
  ylab("Survival rate") + guides(color=guide_legend(title="fatality"), fill=guide_legend(title="fatality")) + 
  coord_flip()
p

# worst airplane
worst_airplane <- casn %>% filter(!is.na(Type)) %>% 
  filter(is_army == FALSE) %>% 
  group_by(Type) %>% 
  summarise(Total_occupants = sum(Total_occupants), Total_fatalities = sum(Total_fatalities),
            Total_survivors = sum(Total_survivors), Survival_rate = 100*Total_survivors/Total_occupants) %>% 
  ungroup() %>% 
  filter(Total_occupants > 500) %>% 
  top_n(20, wt = desc(Survival_rate)) %>% 
  arrange(Survival_rate)

p = ggplot(data = worst_airplane, mapping = aes(x = reorder(Type, Survival_rate), y = Survival_rate, fill = Total_fatalities)) + 
  geom_bar(stat="identity") +
  ggtitle("Worst Airplanes with lowest survival rate") + 
  xlab("Airplane") + 
  ylab("Survival rate") + guides(color=guide_legend(title="fatality"), fill=guide_legend(title="fatality")) + 
  coord_flip()
p

# worst route
worst_departure_airport <- casn %>% filter(!is.na(DepartureAirport)) %>% 
  filter(is_army == FALSE) %>% 
  group_by(DepartureAirport) %>% 
  summarise(Total_occupants = sum(Total_occupants), Total_fatalities = sum(Total_fatalities),
            Total_survivors = sum(Total_survivors), Survival_rate = 100*Total_survivors/Total_occupants) %>% 
  ungroup() %>% 
  filter(Total_occupants > 500) %>% 
  top_n(20, wt = desc(Survival_rate)) %>% 
  arrange(Survival_rate)

p = ggplot(data = worst_departure_airport, mapping = aes(x = reorder(DepartureAirport, Survival_rate), y = Survival_rate, fill = Total_fatalities)) + 
  geom_bar(stat="identity") + scale_fill_gradient(low="midnightblue", high="darkred") +
  ggtitle("Worst Departure Airports with lowest survival rate") + 
  xlab("Departure Airport") + 
  ylab("Survival rate") + guides(color=guide_legend(title="fatality"), fill=guide_legend(title="fatality")) + 
  coord_flip()
p


مدل اولیه دسته بندی دلایل سقوط پروازها
برای این منظور می خواهیم از مدل lda استفاده کنیم. به همین دلیل ابتدا از داده های اصلی تنها شماره ی تصادف و Narrative را انتخاب می کنیم. سپس علت سقوط را به کلمات آن تبدیل می کنیم و stopwords را از آن حذف می کنیم. در نهایت نیز تعداد تکرار هر لغت را برای هر Narrative بدست می آوریم. سپس از آنجایی که LDA با DocumentTermMatrix کار می کند، ساختار داده ی خود را به این صورت تغییر می دهیم. در نهایت نیز مدل خود را با ۲۰ topic لرن می کنیم. از آنجایی که لرن مدل وقت گیر است، مدل را برای استفاده ی آینده ذخیره می کنیم.

cause <- casn %>% select(occ_no, Narrative)

# split into words
cause_word <- cause %>%
  unnest_tokens(output = word, input = Narrative)

word_counts <- cause_word %>%
  anti_join(stop_words) %>%
  count(occ_no, word, sort = TRUE)

flight_dtm <- word_counts %>%
  cast_dtm(occ_no, word, n)
accident_lda <- LDA(flight_dtm, k = 20, control = list(seed = 1234))
saveRDS(accident_lda, file="Data/lda.rds")
accident_lda = readRDS(file="Data/lda.rds")
accident_lda
A LDA_VEM topic model with 20 topics.

سپس برای هر topic پنج کلمه ای که بیشترین احتمال حضور در این موضوع دارد را نمایش می دهیم.

accident_topics <- tidy(accident_lda, matrix = "beta")

top_terms <- accident_topics %>%
  group_by(topic) %>%
  top_n(5, beta) %>%
  ungroup() %>%
  arrange(topic, -beta)

top_terms %>%
  mutate(term = reorder(term, beta)) %>%
  ggplot(aes(term, beta, fill = factor(topic))) +
  geom_col(show.legend = FALSE) +
  facet_wrap(~ topic, scales = "free", nrow = 5) +
  coord_flip()

در نمودارها مشاهده می کنیم که گروه های مختلفی از جمله، وضعیت جوی و محیطی، مشکل فنی هواپیما، مشکلات باند فرودگاه، هایجک کردن هواپیما و … قرار دارد. گروه های این مدل دور از انتظار و غیرمنطقی نیست، در نتیجه مشاهده می کنیم که مدل معقولی داریم.(البته در فازهای بعدی روش های بهتری برای ارزیابی مدل به کار می بریم.)

همانطور که در بالا مشاهده می کنیم، کلمات رایجی همچون aircraft و یا crashed نیز در بین کلمات حضور دارند. به همین دلیل احتمالا در فاز بعدی، این مدل را از طریق حذف کلمات مربوط به سقوط بهبود می دهیم.